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ABSTRACT 



We present results of three dimensional simulations 
of the uppermost part of the sun, at 3 stages of its 
evolution. Each model includes physically realistic 
radiative-hydrodynamics (the Eddington approxima- 
tion is used in the optically thin region), varying 
opacities and a realistic equation of state (full treat- 
ment of the ionization of H and He) . In each evolu- 
tion model, we investigate a domain, which starts at 
the top of the photosphere and ends just inside the 
convection zone (about 2400 km in the sun model). 
This includes all of the super-adiabatic layer (SAL). 
Due to the different positions of the three models in 
the log{g) vs logT^ff plane, the more evolved mod- 
els have lower density atmospheres. The reduction 
in density causes the amount of overshoot into the 
radiation layer, to be greater in the more evolved 
models. 



Key words: Radiative hydrodynamics, compressible 
turbulence. 
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Figure 1. Evolutionary track for the sun from the 
pre-main sequence stage to the giant branch 



super adiabatic layer (SAL). 

The SAL is important for the following reasons: 



1. INTRODUCTION 



In most of the convection zone, heat transport can 
be modelled quite well by using the mixing length 
theory. Due to vigorous turbulent mixing, the en- 
tropy is almost constant. This well mixed layer has 
an extremely small super-adiabatic temperature gra- 
dient (V — Vad ~ 10~®) and is practically opaque 
(the optical depth is greater than 10^). In such a 
layer, radiative transport accounts for a very small 
fraction of the heat flow. This changes further up. 
Near the top of the convection zone the gas density 
is very low, so the enthaply flux is not big enough 
to carry all of the outward flowing heat. Radiation 
must carry most of the heat flux. To enable radi- 
ation to transport this heat flux, the temperature 
gradient must increase significantly. This region, in 
which V ~ Vad is of order unity, is known as the 



• It is the site of maximum amplitude of the 
p-mode oscillations (particularly the high fre- 
quency modes). 

• As most of the entropy change occurs within this 
layer, it plays an important role in determin- 
ing the radius of stellar models (Larson 1974). 
(Note dS/dlnP = Cp[V - Vad]-) 

• The amount of overshoot into the enveloping 
photosphere depends on the structure of the 
SAL. 



In this poster we will describe three simulations of 
stellar atmospheres each centered around the SAL. 
We consider the sun at three different stages of its 
evolution, namely the ZAMS, present sun and sub- 
giant. These are shown as three points in the theo- 
retical HR diagram, figure |l|. 
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2. MATHEMATICAL MODEL 



2.1. Hydro dynamic equations 



In the outer layers of the sun the Mach number, 
{v/vg)"^ can be of order unity (Cox and Giuh 1968) 
(v is the flow velocity and Vg is the sound speed) . In 
such an environment, the governing hydrodynamic 
equations are the fully compressible Navier Stokes 
equations (cf, Kim et al 1995). These are modified 
to include radiative energy transport by the insertion 
of the Qrad term (described later) into the energy 
equation. The full set of governing equations are: 
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dpv/dt 
dE/dt 
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where E = e + pv^/2 is the total energy density 

and p, V, P, e and g, arc the density, velocity, pres- 
sure, specific internal energy and acceleration due 
to gravity, respectively. Ignoring the coefhcicnt of 
bulk viscosity (Becker 1968), the viscous stress ten- 
sor for a Newtonian fluid is Sjj = ^{dvi/dxj + 
dvj/dxi) — 2/i/3(V • 'v)5ij. In LES, the dynamic vis- 
cosity /i is increased so that it represents the effects 
of Reynolds stresses on the unresolved or sub-grid 
scales (Smagorinsky 1963), 

^i = p{c^^f{2a^.aY/\ 

The colon inside the brackets denotes tensor con- 
traction of the rate of strain tensor (Jij = (ViVj + 
W jVi)/2. The SGS eddy coefficient c^, is set to 0.2, 
the value for incompressible turbulence, and A is an 
estimate of the local mesh size. The present formula- 
tion ensures that the grid Reynolds number Axv/v 
is of order unity everywhere. To handle shocks, fi 
is multiplied by 1 + C • (V ■ v)^, where the constant 
C is made as small as possible to maintain numer- 
ical stability. As fx is dependent on the horizontal 
divergence, any large horizontal velocity gradients 
are smoothed out by the increased viscosity. The 
diffusive flux f = — (/i/Pr)rVS', where the horizon- 
tal mean of WS < i.e. the convection zone, and 
f = — (/xCp/Pr)VT where the horizontal mean of 
V5 > i.e. the radiation zone. The Prandtl number 
Pt = iz/k {v is the kinematic viscosity and k, is the 
thermal diffusivity) is 1/3. Due the inclusion of ra- 
diative energy transport, the effective Pr is actually 
much smaller and not constant. 



2.2. Radiative energy transport 

In the deeper part of the domain, radiative transfer 
is treated by the diffusion approximation. 



where k, is the Rosseland mean opacity, a is the Boltz- 
mann constant and c is the speed of light. 

In the shallow region, the photon mean free path is 
at least one tenth of the depth of the atmosphere so 
the diffusion approximation may not apply. Instead 
<3rad is computed as 

<9rad = 4kp( J - B) 

where mean intensity J is computed by using the 
generalised three-dimensional Eddington approxima- 
tion (Unno and Spiegel 1966), 
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V J - KpJ + npB = 0. 



Q 



rad 



3kp 



VT 



This is exact for isotropic radiation, and in non- 
equilibrium the Eddington approximation describes 
the optically thick and thin regions exactly. 



3. NUMERICAL METHODS 



The simulation domain is a small box of aspect ratio 
1.5, which includes the photosphere and top part of 
convection zone. In the sun model this represents 
6.936 X 10^°cm < R< 6.960 x 10^°cm. Each of the 
three models span about 8 pressure scale heights. 

The governing equations are discretiscd in cartesian 
coordinates on 80^ uniformly spaced grid points. Us- 
ing a code developed by Chan and Wolff (1982), an 
implicit scheme (the Alternating Direction Implicit 
Method on a Staggered grid or ADISM) relaxes the 
fluid to a self consistent thermal equilibrium. In the 
fully relaxed layer, the energy flux leaving the top 
of the box is within 5 % of the input flux at the 
base. Next a second order explicit method (Adams 
Bashforth time integration) gathers the statistics of 
the time averaged state. The statistical integration 
time is about 2500 seconds of solar surface convec- 
tion, and requires about a million time steps. On an 
667 MHz Alpha processor, each integration step re- 
quires about 10 seconds of CPU time. Consequently 
each simulation takes at least 3 months to run. 



3.1. Modelling a stellar atmosphere 

A Standard Solar model, calculated with the Yale 
Stellar Evolution model (Guenther & Demarque 
1997), is used to compute the initial stratification 
i.e. run of pressure, temperature, density, internal 
energy, for the box. The model uses realistic physics, 
Alexander low temperature opacities and the OPAL 
opacities and Equation of State (Hydrogen and he- 
liTim ionisation zones are included). The hydrody- 
nainical simulations use identical opacities and e.o.s, 
as used in the 1-d stellar model (Kim & Chan 1998). 
The horizontal boundaries are periodic, while the 
vertical boundaries are stress free. A constant heat 
flux flows through the base, and the top is a perfect 
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conductor. To ensure mass, momentum and energy 
are fully conserved, we use impenetrable (closed) top 
and bottom boundaries. 



4. RESULTS AND INFERENCES 



4.1. Turbulent quantities 

Each parameter, which consists of a mean and a fluc- 
tuating part, is computed from the 3-d statistically 
averaged flow. For a given parameter X, the turbu- 
lent part is approximated by the variance, 

where the overbar denotes horizontal and temporal 
averaging, and X is the total quantity (mean plus 
fluctuating). The autocorrelation between two fluc- 
tuating quantities X[ and is computed as 



C[X[X!,] = (Xi - Xi)((X2 - X2)/xtX2. 



4.2. Structure of the SAL in the three models 

Figures 2, 3 and 4 show the entropy gradient and the 
autocorrelation between vertical velocity fluctuation 
and temperature, C[V^'T'], vs depth, using InP as 
the unit of depth. 

The SAL peak is higher and sharper in the subgiant 
model compared to the ZAMS model. This is be- 
cause density is much lower in the SAL region for the 
subgiant, than it is for the ZAMS. Convection is less 
efficient, and thus V — Vad (or similarly dS/dlnP) 
needs to be steeper and higher, so that radiation can 
carry more heat. Generally, the SAL regions are ra- 
dially further outwards in progressively more evolved 
models. 

The difference between the top of the convection 
zone, as defined by dS/dz = (Schwarzschild 
criterion), and the point where the enthapy flux, 
pc-pC\V!,T']'wt^ is zero, is one measure of the amount 
of overshoot into the photosphere. The figures 
(which exclude the two outermost levels) show that 
the overshoot increases slightly between the ZAMS 
and sun models, while extending high up into the 
lower density region in the subgiant model. It is 
important to note that bcx;ause tlic^ grid points are 
uniformly spaced, the resolution of the SAL is worse 
in the subgiant than in the ZAMS. As the SAL oc- 
curs higher up in the subgiant where the pressure 
scale height is smaller, there are less grid points be- 
tween InP equals 11 and 12, than between 12 and 
13. Clearly it is very important to have enough grid 
points to model the SAL in more evolved stellar mod- 
els. 
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Figure 2. Horizontally averaged entropy gradient and 
correlation between vertical velocity and temperature 
vs depth, for the ZAMS model. The crosses represent 
the actual grid points. The normalised correlation is 
rescaled to fill the same figure. 
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Figure 3. Same as figure 2 for the sun model. 
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Figure 4- Same as figure 2 for the subgiant model. 
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4.3. Turbulent velocities 
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Figure 5. For the ZAMS model, the turbulent ve- 
locities in the horizontal (u and v) and vertical (w) 
directions, versus depth. The closeness of the two 
horizontal velocities confirms that the simulation is 
statistically converged. 

Present Sun 



3xlC 



2x10" 



1 1 1 1 1 1 1 


1 1 1 1 _ 


^ %^ 


u(diamonds) I 
v(crosses) I 
w(pluses)^ ~ 






- + 




_ + ^ 




I + 


+ z 


- + 


+ z 


~ 1 1 1 1 1 1 1 





14 

In P 



18 



Figure 6. Same as figure 5 for the sun model. 
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The turbulent vel ociti es u, v and w, computed as de- 
scribed in section 4.1 are shown in figures 5,6 and 7. 
The closeness of the turbulent horizontal velocities 
(m and v) indicates the degree of statistical conver- 
gence. The ZAMS and sub-giant models are fully 
converged, while the sun is close to convergence. 
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Figure 7. Same as figure 5 for the subgiant model. 



